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Abstract. The OPERA detector at the Gran Sasso underground laboratory (LNGS) was used to mea- 
sure the atmospheric muon charge ratio i?^ = /N^- in the TeV energy region. We analyzed 403069 
atmospheric muons corresponding to 113.4 days of livetime during the 2008 CNGS run. We computed 

separately the muon charge ratio for single and for multiple muon events in order to select different energy 
regions of the primary cosmic ray spectrum and to test the dependence on the primary composition. 
The measured Rfi values were corrected taking into account the charge-misidentification errors. Data have 
also been grouped in five bins of the "vertical surface energy" f ^ cos 6. A fit to a simplified model of muon 
production in the atmosphere allowed the determination of the pion and kaon charge ratios weighted by 
the cosmic ray energy spectrum. 
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1 Introduction 

Primary cosmic rays (typically protons) impinging on the 
Earth's atmosphere produce showers of secondary parti- 
cles which propagate down to the ground level. Most of 
the interaction products are tt and K mesons which in 
turn decay or interact, depending on their energy and on 
the air density profile they pass through. The decay of 
tt" mesons gives rise to the electromagnetic component 
of the showers, the decay of tt^ and mesons yields 
mostly muons which are the most penetrating charged 
particles and therefore the most abundant charged com- 
ponent at sea level. In particular only the most energetic 
muons can penetrate deep underground. The Gran Sasso 
laboratory (LNGS) is located at an average depth of 3800 
m.w.e. and the minimum muon energy required to reach 
the underground depth is around 1.5 TeV while the resid- 
ual underground energy is about 270 GeV averaged over 
all directions [1 . 

The muon charge ratio = A'p+ /N^- defined as 
the number of positive over negative charged muons, re- 
sults from several contributions; the primary cosmic ray 
composition (in particular the ratio of protons over heav- 
ier primaries), hadronic-interaction features, atmospheric 
conditions (negligible above a few GeV) and, at very high 
energy, the contribution of muons from charmed particle 
decays (prompt muons) 2 . The muon charge ratio at sea 
level was extensively studied in the past since it is an in- 
dicator of important aspects of cosmic rays and particle 
physics. 



An exhaustive compilation of measurements in a wide 
energy range is reported in Ref. . In the interval from a 
few hundred MeV to 300 GeV the muon charge ratio 
stays around 1.27. At higher energies several competing 
processes can affect its value. Since strong interaction pro- 
duction channels lead to a / K" ratio higher than the 
tt^/tt" one and the fraction of muons from kaon decays in- 
creases with the energy, the muon charge ratio is expected 
to rise as the energy increases. On the other hand, as the 
zenith angle increases and hence longer lived mesons have 
a higher probability to decay in the less dense layers of 
the high atmosphere the fraction of muons from pion de- 
cay increases and the muon charge ratio decreases. We 
also expect a dependence of the muon charge ratio on the 
underground muon multiplicity , which is related to the 
energy of the primary cosmic rays and to their chemical 
composition. For primaries different from protons the pos- 
itive charge excess is reduced and so is the muon charge 
ratio [4]. 

A simplified model of the atmospheric muon charge 
ratio is obtained from the muon spectrum [S] 

Npar 

^ 1-Zjvw l + h£Je,{e) ^ ' 

where 'pN{£fj.) — ^o^ii is the primary spectrum of 

nucleons (evaluated at the muon energy in the atmosphere 
£fi) with a spectral index 7 + 1 ~ 2.7. For each of the Npar 
muon parents (tt, K, charmed particles etc.) the constants 
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ai and bi contain the kinematical factors for the decay 
into muons, ei{9) are the critical energies defined as the 
energies above which interaction processes dominate over 
decay. They depend on the ratio rrii / Ti (mass over rest life- 
time of the muon parent) and on the atmospheric profile 
density and therefore on the zenith angle 9. A good ap- 
proximation for ei{9) which takes into account the Earth's 
curvature is 

- ^ (2) 



cos 9* 



with 



cost 



Re 



Re + h 



(3) 



where Re is the Earth's radius and h is the muon produc- 
tion height. By using eq.[3]the zenith angle is evaluated at 
the muon production point and not at the detector site. 
By choosing h = 30 km an agreement within 5% with the 
precise £i{9) computation is obtained 6 . 

The spectrum weighted moments are defined as: 



Z,.. 



1 daij 

CTi, dxial 



■{xiabT ^dxiab 



(4) 



where aij is the inclusive cross-section for the production 
of a particle j from the collision of a particle i with a 
nucleus in the atmosphere and xiab — Ej / Ei is the energy 
fraction carried by the secondary particle. Use of the Z 
factors shows explicitly that particle production in cosmic 
ray cascades is concentrated in the forward fragmentation 
region at large xiab- 

From Eq. [T] the i-th contribution to the muon flux is 
proportional to Z^i and is suppressed for :s> ei{9). 
Therefore each contribution to the muon charge ratio pro- 
duced by different muon parents can be disentangled by 
studying the muon charge ratio as a function of the muon 
energy 5^. 

Considering only tt and K parent mesons Eq. [TJcan be 
separated for and /i^ 



*P ± OC 'T' 

^ \l + bT,£^cos9*/e^ 1 + bK£^cos9* /ex 



(5) 



where e^r — 115 GeV and ex — 850 GeV are the pion 
and kaon critical energies along the vertical direction, re- 
spectively. These energies have to be compared with the 
corresponding value for charmed particles, ex > 10^ GeV. 
The prompt muon component (from charmed particles) is 
therefore isotropically distributed since the corresponding 
cos 9* factor is suppressed, at least in the TeV region. 

Eq. [5] contains most of the aspects already discussed. 
First we note that the correct variable to describe the 
evolution of the charge ratio is the product f ^ cos 9* , the 
"vertical surface energy" [71IH]. Moreover the two energy 
scales which determine the pion and kaon contributions 
to Rf^ are the critical energies and e^- The evaluation 
of the muon surface energy £^ depends on the rock depth 
crossed by the muon to reach the detector and therefore 
the distribution of £^ cos 9* is related to the shape of the 
overburden. Measurements of the muon charge ratio at 



high energies and large zenith angles, corresponding to 
{£f_iC0s9*) ~ 0.5 TeV, are given in Ref. |9]. More recent 
data with large statistics at cos 9*) ~ 1 TeV are pre- 
sented in Ref. [TU]. These results suggest a smooth tran- 
sition toward the energy region where kaon contribution 
becomes significant. 

The LNGS laboratory is located at (£^cos0*) ~ 2 
TeV, well above the kaon critical energy ex- This allows 
the measurement of the ratio Z]^j^+ jZj^j^- whose value is 
poorly known in the fragmentation region. This has also a 
strong impact, for instance on the evaluation of the flux of 
TeV atmospheric neutrinos, which are dominated by kaon 
production. OPERA is the first large magnetized detec- 
tor that can measure the muon charge ratio at the LNGS 
depth, with an acceptance for cosmic ray muons coming 
from above A — 599 m^- sr {A = 197 m^ sr for muons 
crossing the spectrometer sections). 

The paper is organized as follows. In Sec.[2]the detector 
is briefly described, while in Sec. [3] we describe data selec- 
tion and reconstruction, Monte Carlo simulation and data 
reduction. In Sec. H] the muon charge ratio at the LNGS 
underground depth and its systematic error are evaluated. 
Finally in Sec. [5] we give R^ as a function of the under- 
ground momentum and of the variable cos 9* fitted to 
Eq.El 



2 The OPERA detector 

OPERA is a hybrid experiment with electronic detectors 
and nuclear emulsions located in Hall C of the under- 
ground Gran Sasso Laboratory in central Italy [TT]. The 
main physics goal of the experiment is to observe neutrino 
flavor oscillations through the appearance of neutrinos 
in the CNGS beam [T^l. The detector design was op- 
timized to identify the r lepton via the topological obser- 
vation of its decay: this requires a target mass of more 
than a kton to maximize the neutrino interaction proba- 
bility and a micrometric resolution to detect the r decay. 
To accomplish these requirements the detector concept is 
based on the Emulsion Cloud Chamber (ECC) technique 
combined with real-time detection techniques ( "electronic 
detectors" ) . 

The ECC basic unit in OPERA is a "brick" made of 
56 lead plates, 1 mm thick, providing the necessary mass 
to cope with the small neutrino cross-section, interleaved 
with 57 nuclear emulsion films (industrially produced), 
providing the necessary spatial and angular resolution to 
identify tau decay topologies. In total, 150000 bricks have 
been assembled reaching the overall mass of 1.25 kton. 

The electronic detectors are used to trigger the neu- 
trino interactions, to locate the brick in which the interac- 
tion occurred, to identify muons and measure their charge 
and momentum. 

The detector is composed of two identical parts called 
supermodules (SM), each one consisting of a target sec- 
tion and a magnetic spectrometer. In the target the bricks 
are arranged in 29 vertical "walls" transverse to the beam 
direction interleaved with electronic Target Tracker (TT) 
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Fig. 1: Schematic view of a charged particle crossing one spectrometer. The six PT stations are shown in dark grey; the 
24 iron slabs (12 per arm) interleaved with 22 RPC planes are shown in light grey. Each spectrometer arm provides an 
independent measurement of charge/momentum, provided the track is reconstructed in at least one station (or station 
doublet) in each side of the arm. 



walls. Each TT wall consists of a double layer of 256 scin- 
tillator strips, for vertical and horizontal coordinate mea- 
surements. The TTs can trigger the data acquisition and 
locate the brick in which the interaction occurred. 

The target section is followed by a magnetic spectrom- 
eter (see Fig. [1]). A large dipolar iron magnet is instru- 
mented with Resistive Plate Chambers (RPCs). The mag- 
netic field is 1.53 T, directed vertically transverse to the 
neutrino beam axis. The RPC planes are inserted between 
the iron slabs: they provide the tracking inside the magnet 
and the range measurement for stopping muons |13) . 

The deflection of charged particles in the magnet is 
measured by six stations of vertical drift tubes, the Preci- 
sion Trackers (PT), grouped in 3 pairs placed upstream of 
the first, in between and downstream of the second magnet 
arms (Fig. [Ij . Each PT station is formed by four staggered 
layers of aluminum tubes, 8 m long, with 38 mm outer di- 
ameter. The spatial resolution is better than 300 /um in 
the bending (horizontal) plane: this allows the determi- 
nation of the muon charge sign with high accuracy, and 
the momentum measurement with a resolution of better 
than 20% for momenta < 50 GeV/c for charged particles 
coming from the CNGS direction [T3]. The PT system is 
triggered by the RPC timing boards with a configuration 
optimised to collect both beam and cosmic ray muons with 
high efficiency. 

In order to remove ambiguities in the reconstruction 
of particle trajectories, in particular in multi-track events, 
each spectrometer is instrumented with additional RPC 
planes (XPC), with two crossed strip planes rotated by 
±42.6° with respect to the horizontal. 

Two RPC planes (VETO) are placed in front of the 
detector, acting as a veto for charged particles originating 
from the upstream material (mainly muons from neutrino 
interactions in the rock). 



Finally we emphasize that for this analysis the OPERA 
detector was used differently from what it was conceived 
for. This is particularly true for the PT system which was 
configured and optimized to reconstruct and measure par- 
ticles traveling along the CNGS direction. 



3 Data analysis 

3.1 Data collection and selection 

The results presented are based on data recorded during 
the CNGS physics run, from June 18th until November 
10th 2008. The detector ran in the standard configura- 
tion, with the magnetic field directed along the vertical 
axis in the first arm of both spectrometers, and in the op- 
posite direction in the second arm. Moreover a sample of 
cosmic ray muons was collected with the magnetic field 
switched off in order to improve the alignment between 
PT stations and to evaluate systematic uncertainties. A 
limited data sample was obtained inverting the magnet 
polarity to cross check the charge reconstruction. 

The data acquisition was segmented in "extraction pe- 
riods" of about 12 hours each. Only periods of data tak- 
ing where all the main detector subsystems ran in stable 
conditions were considered. They amount to about 78% of 
the total duration of the run. The total number of selected 
events is 403069 corresponding to 113.4 days of livetime. 

3.2 Event reconstruction 

The OPERA standard software for beam event recon- 
struction was complemented with a set of dedicated soft- 
ware tools developed for cosmic ray events. Once the event 
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is tagged as "off-beam", that is outside the CNGS spih 
window it is classified as cosmic and processed in a dedi- 
cated way. This choice was required by the different topolo- 
gies of beam and cosmic ray events. Beam events come 
from a well defined direction (the CNGS one) and the re- 
construction code is optimized to follow a single long track 
(the muon escaping from the neutrino-interaction region) 
on the z-axis. Cosmic ray events come from all directions, 
they are not generated within the target and a fraction 
of them (~ 5% in OPERA) are muon bundles. A brief 
description of the code follows. 

The reference frame is defined to have the z axis along 
the Hall C longitudinal direction (from north to south), 
y perpendicular to the floor pointing toward the zenith 
and X describing a right-handed frame. In this coordinate 
system, the zenith direction 6 is defined by the angle with 
the y axis, the azimuth direction (jj by the angle with the z 
axis. Event reconstruction is performed separately in the 
two projected views Txz and Tyz. First the event direc- 
tion is determined by using the Hough transform. Using 
a Monte Carlo simulation we estimated an angular reso- 
lution of better than 0.5°, both in the 9 and 4> directions, 
for single as well as for multiple track events. Then, the 
direction information is used to subdivide the T^z and Tyz 
views in slices 25 cm wide having the same slant as the 
reconstructed direction. The hits within the same slice are 
then processed separately to search for a track "seed" of 
at least three aligned points. If a seed is found, all the 
other hits in the corresponding projected view are linked 
to the selected track according to pre-defined tolerances. 
Tracks independently reconstructed in each view are then 
merged together to build the three-dimensional event. 

For this analysis, particular attention was devoted to 
the reconstruction of tracks with the Precision Tracker. 
A description of the PT system is available in [H], while 
the reconstruction procedures are detailed in T5| . Here we 
briefly mention the main steps used to extract charge and 
momentum from PT hits. 

A muon crossing the spectrometer is deflected in the 
horizontal plane (Fig. [1]) . Let be the angle between the 
particle direction and the z axis, then the deflection Acj) 
is the difference between the two angles measured at both 
sides of each magnet arm. Each (j) value is obtained by 
fltting the spatial information provided by the PT system 
which is made of 8 layers of drift tubes arranged in two 
"stations" . Since cosmic ray tracks are almost uniformly 
distributed in (f), many of them traverse only one single 
statiorQ of a pair and we refer to them as singlets, other- 
wise we call them doublets. 

Due to the different lever arm the angular resolution 
for singlets is worse than for doublets. For tracks par- 
allel to the z-axis {(f> = 0) it is cr^ ~ 1 mrad for sin- 
glets and ~ 0.15 mrad for doublets. Since for tilted tracks 
the number of fired tubes and their mutual distances are 
larger than for tracks with cj) = 0, the errors on the slopes 
decrease. In order to increase the statistics we decided 
to consider both singlets and doublets: the percentage of 



cases in which both angles are reconstructed from doublets 
is ^55% of the total, ~9% are from singlets and the re- 
maining 36% are from mixed configurations, namely cases 
where an angle is reconstructed from a doublet and the 
other angle from a singlet. 

Considering tracks with = and the Multiple Coulomb 
Scattering (MCS) within one magnet arm the total uncer- 
tainty on A(j) is 



o.oisey _d_ 
) x'o 



p 



(6) 



where d = 0.6 m is the magnetized iron thickness, Xq = 
0.0176 m is the iron radiation length and p is expressed in 
GeV/c. Since the defiection due to the magnetic field is 



O.Wd 

P 



(7) 



where B = 1.53 T, the requirement A(j)/<T/^^ > 1 pro- 
vides an estimate for the maximum detectable momen- 
tum, pmax — 1-25 TeV for doublets, Pmax — 190 GeV for 
singlets and Pmax — 260 GeV for mixed configurations. 

For muon momenta p <C Pmax the measurement er- 
ror can be neglected and the only contributions to the 
Z\0 uncertainty come from the MCS. In this ideal case 
the ratio A4>b / A(j)Mcs 3.5 corresponds to a charge- 
misidentification ij (defined as the fraction of tracks re- 
constructed with wrong charge sign) below 10"'^. In real- 
ity there are other effects which spoil the resolution and 
therefore the charge identification capability, as detailed 
in Sec. O 

To measure the charge and the momentum of a particle 
at least one Acj) angle is needed (for tracks parallel to the 
z-axis there can be up to 4 independent angles). For each 
reconstructed A(j)i the track momentum is computed by 
using the formula |15) 



Pz = 



l{dE/dz) 



1 - eyL-p[A(l),{dE I dz) / eB] 




(8) 



^ A PT station is made of 4 layers of drift tubes. A track is 
reconstructed when Ntubes > 4. 



where I — 0.82 m is the total arm length (including RPC 
gaps), B = Bd/l is the effective magnetic field and dE/dz 
is the ionization energy loss in the magnetized iron (which 
depends logarithmically on the muon momentum). The 
term in the square root takes into account the track slope 
Syz in the Tyz view. The muon charge is determined from 
the sign of the Z\0i angle, accounting for the particle ar- 
rival direction and the field orientation in the arm. The 
final muon momentum and charge are computed as the 
weighted average of the independent measurements. 



3.3 Monte Carlo simulation 

Two different Monte Carlo simulations have been used in 
order to meet two different and conflicting requirements. 
From one side, one needs a large statistical sample of un- 
derground cosmic ray muons for calibration purposes and 
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angle (rad) N' 

(a) Functional form of M{(l)) (b) Rescaled data 

Fig. 2: Cut on the number of fired PT tubes/station. The geometrical dependence of the number of fired tubes on the 
(/> angle is shown in (a), together with a polynomial fit; the rescaled distributions are shown in (b). A 3cr cut of the 
Gaussian fit to Monte Carlo events where secondary particle production was switched off, was applied to the rescaled 
data (see text). 



to correct for detector effects. This task was accomplished 
with a code which generates muons directly at the detec- 
tor level fast enough to produce the required statistical 
amount of events. On the other side, information on cos- 
mic ray primaries and the links between underground and 
surface variables required the use of a complete generator 
which simulates the full cosmic ray cascade in the atmo- 
sphere. 

The first Monte Carlo code (MCI hereafter) is based 
on a fast parameterized event generator developed in the 
framework of the MACRO experiment [TBI and adapted 
for OPERA. The generator considers the MACRO pri- 
mary cosmic ray composition model |17j and, for each 
primary mass, it extracts from built-in probability tables 
the underground muon multiplicity. This choice ensures a 
self-consistency on the predicted muon fiux underground 
since the composition model and the probability tables 
were obtained using the same hadronic interaction model. 
This means that the systematic errors on the primary 
composition and on the interaction model cancel and the 
Monte Carlo generator predicts the correct muon flux in 
the Hall B of LNGS. The event direction {9, (f) is sam- 
pled, together with the muon radial distribution with re- 
spect to the shower axis and the muon underground mo- 
mentum. Finally, the event kinematics is processed within 
the OPERA standard software chain, from the generation 
of the relevant physical processes in the electronic detec- 
tors up to the reconstruction level. The atmospheric muon 
charge ratio is introduced by hand {R^ = 1.4). Taking 
into account the livetime normalization, the ratio between 
OPERA data and Monte Carlo MCI predictions is 



Rate_R_EAL 
RatcMci 



(95.9 ±0.3)% 



(9) 



The difference from unity of the ratio is mainly due to 
the subdetectors efficiency mismatch between data and 
detector simulation. 

The second Monte Carlo program (MC2) is based on 
the more detailed and physics-inspired Monte Carlo pro- 
gram described in Ref. [TH], based on the FLUKA code 
[TO] . Here all the main physical processes are implemented, 
from the primary cosmic ray interactions and shower prop- 
agation in the atmosphere up to the muon transport in the 
overburden. The primary composition model is described 
in Ref. [50] based on a global fit of several experimental 
observations. This code is predictive of the muon charge 
ratio and allows the study of the relations between under- 
ground and surface muon energies. Due to its complexity 
the event production with MC2 is statistically limited. 

In the following with the term "Monte Carlo" we refer 
to MCI unless otherwise specified. 



3.4 Data reduction 

A set of data quality cuts were made in order to isolate a 
clean sample of reconstructed muon events. First at least 
one reconstructed A(j) angle is required for each event {ac- 
ceptance cut). Then events with a large number of PT 
hits potentially dangerous for the muon charge determi- 
nation are removed [clean PT cut). This typically occurs 
when some drift tubes are fired by secondary particles 
(5-rays, showers etc.) and the best track could result 
from a fake tube configuration. In order to evaluate the 
maximum number of fired tubes/track allowed by geomet- 
rical considerations a special version of MCI switching off 
delta ray and secondary particle production was run. By 
naming M and N the number of fired tubes from Monte 
Carlo simulation and experimental data, respectively, we 
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derived the functional form M = M{(f)), a six-order poly- 
nomial shown in Fig. [2al M{(j)) was used to rescale the 
experimental distribution N as N = N — M{(j)) (Fig. l2b|) . 
We considered only tracks with N < 3(j (one-sided cut), 
where a is the standard deviation of the Gaussian fit to 
N . We verified by visual inspection that events rejected 
by the latter cut are characterized by a large number of 
additional fired tubes in the neighbourhood of the correct 
ones. 

A further cut was applied on the Acj) angle {deflection 
cut). Events having a A(j) smaller or compatible with the 
experimental resolution were rejected. On the basis of the 
plot shown in Fig. [3] where as expected for small deflec- 
tion values Rf^ — 1, only events with A(j)/aA4> «^ were 
selected. The effect of this cut is visible in Fig. 2] in which 
the Acj) distribution is shown before and after its applica- 
tion. In these plots, experimental data (black points) are 
plotted with the corresponding Monte Carlo distributions 
split in the two regions corresponding to positive particles 
(Qtrue > 0) and negative particles {qtme < 0) . The charge- 
misidentification rj corresponds to the overlapping region 
of the two distributions. Averaged over all the event sam- 
ples T] is reduced from 0.080±0.002 to 0.030±0.001 by this 
cut. The source of events with large A<j> angles and recon- 
structed with wrong charge-sign was investigated. A visual 
scan of Monte Carlo events confirms the hypothesis that 
they are due to secondary particles in the neighbourhood 
of the true muon track: if the two tracks are very close, it 
may happen that the track reconstructed with the best 
is the wrong one. A further selection Acj) < 100 mrad was 
used to reject these fake tracks, with a small impact on 
the statistics. This last selection affects the sample with 

< 5 GeV/c. 
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Fig. 3: Dependence of the measured charge ratio on the 
deflection angle expressed in units of experimental reso- 
lutions. A cut at A(j)/aA<p > 3 was applied in the data 
analysis. Note that the fitted value i?^ = 1.342 ± 0.015 
was obtained with the bins indicated in the plot (the first 
3 bins have not been used). 



The muon charge ratio was computed separately for 
single muon events (i.e. event multiplicity = 1) and 
multiple muon events (n^ > 1). Single muon events are 
selected by requiring single tracks in each projected view 
merged in the three dimensional space. Multiple muon 
events are selected by requiring a muon multiplicity > 2 
in both views, with tracks identified and unambiguously 
merged in 3D space. 

Tab. [1] lists the number of events remaining at each 
stage of the selection process. Note that data and MCI 
event rates are absolute (given in day~^) and not normal- 
ized one to the other. Also note that the clean PT cut has 
a stronger impact on data reduction and that the effect 
on the experimental data is different from that of Monte 
Carlo. This was expected since the percentage of events 
with PT digits not related to the muon track is intrinsi- 
cally larger in the experimental data. The clean PT cut 
was tuned in order to be left with a clean data sample at 
the expense of a considerable loss of statistics. 

3.5 Alignment of the PT system 

The measurement of the muon charge is strongly affected 
by the alignment precision of the PT system. Misalign- 
ment effects have "global" or "local" contributions. To 
correct for global effects, which are the dominant ones, 
each station is treated as an independent rigid body and 
relative rotations and translations of one station with re- 
spect to the others are searched for. The local misalign- 
ment contribution takes into account possible distortions 
or bendings within each station. 

A first alignment campaign was carried out with a 
theodolite to measure the position of the PT walls in 
the OPERA coordinate system. Recently a more refined 
alignment using cosmic ray muons was performed. The 
alignment procedure was carried out in two steps a) PT 
stations forming a doublet were aligned with the whole 
data sample, since the space in between has no magnetic 
field and tracks do not suffer any deflection; b) each dou- 
blet (pair) treated as a unit, separated by the iron magnet 
arm, was aligned using special runs with the magnetic field 
switched off. This procedure allowed aligning two PT sta- 
tions within a doublet with a spatial accuracy of ~0.1 mm 
and an angular accuracy of ~0.1 mrad, and to align two 
doublets with an angular accuracy of 0.2 mrad. 

Local effects, such as bendings or distortions, contribute 
at the second order level and due to the present limited 
statistics have not been corrected for. However in Sec. 14.21 
we provide an estimate of the systematic uncertainty on 
i?p introduced by these effects. 

4 Underground muon charge ratio 
4.1 Computation of i?^ 

R^ was computed separately for single and multiple muon 
events. Tab. [2] refers to single muon events where the num- 
ber of positive and negative muons, their ratio, the charge- 
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A({> (mrad) A<|) (mrad) 

(a) Before deflection cut (b) After deflection cut 

Fig. 4: Effect of the deflection cut on A<j) distributions reconstructed exclusively from doublets (shown in the range 
-10 ~ 10 mrad). On the left is shown the distribution before the cut, where the two peaks corresponding to /i^ and 
/z~ are already clearly visible. Black points correspond to experimental data, hatched histograms to Monte Carlo 
simulations, split in the two components qtme > and qtme < 0. The same distributions are shown on the right after 
the application of the deflection cut A(j)/aA4> 3. The overlapping region of the two hatched histograms corresponds 
to the charge-misidentified tracks. 







Data 






MCI 






evt/day 


fi 


f2 


evt / day 


fi 


f2 


Acceptance 


992 


100.0% 




1222 


100.0% 




Clean PT 


515 


51.9% 




959 


78.5% 




Deflection 


391 


39.4% 


76.0% 


708 


58.0% 


73.8% 


Single /i 


379 


38.2% 


96.9% 


673 


55.1% 


95.1% 


Multiple fi 


12 


1.2% 


3.1% 


35 


2.9% 


4.9% 



Table 1: Progressive reduction of the number of events per day after each selection cut, for data (left) and for MCI 
(right). The effect of data reduction is also shown by reporting the fraction of events referred to the original sample 
(fi) and to the previous cut (f2). The total number of experimental events surviving the cuts is 44370. 



misidentification rj and the unfolded charge ratio are re- 
ported. The T] value, defined as the fraction of tracks re- 
constructed with wrong charge sign, was extracted from 
the Monte Carlo simulation. Once rj is known, the un- 
folded charge ratio is obtained according to the formula 
(see Appendix IX)) : 



K 



unf _ 



(1 - r])BJ^'"'' - 77 

_^^meQS + (1 _ 



(10) 



The single muon sample was subdivided into three classes: 
tracks reconstructed exclusively as doublets, tracks recon- 
structed exclusively as singlets and as mixed. We verified 
that the fraction of these classes for experimental data and 
for Monte Carlo simulation are compatible: 54.8% (dou- 
blets), 9.0% (singlets) and 36.2% (mixed) for real data to 
be compared to 52.5%, 10.0% and 37.5% for Monte Carlo 
simulation (the errors are < 0.5%). The final charge ra- 
tio value for single muon events, integrated over all the 
classes, is: 



■Rr^K 1) = 1.377 ±0.014 



(11) 



The same procedure was applied to multiple muon events. 
We selected events with > 1 and reconstructed the 
charge of muons crossing the spectrometer section. Events 
were classified in this category provided that more than 
one muon was reconstructed in the detector even though 
only one charge was measured. In other words, the muon 
multiplicity is used to "tag" events generated by heavier 
and more energetic primaries. The possibility to compute 
the muon charge ratio within the same event is presently 
excluded by the lack of statistics of high multiplicity muon 
bundles. 

The charge ratio is i?™'=°''(n^ > 1) = 919/753 = 1.22 ± 
0.06 and the corresponding unfolded value, obtained from 
Eq.fTU] 

('V > 1) 1-23 ± 0.06. (12) 

This value is 2.4cr away from the value for single muon 
events, consistent with the hypothesis of dilution of 
due to the neutron enhancement in the primary nuclei. 

Tab. [3] gives information obtained with MC2 on some 
variables of single muon events and muon bundles in the 
OPERA detector. In particular, are given the average pri- 
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mary mass number (A) , the average primary energy / nucleo: 
(E/A), the fraction of Hydrogen nuclei over the total (H 
fraction) , the ratio of protons over neutrons in the primary 
radiation Np/Nn and finally the measured muon charge 
ratio 

4.2 Systematic uncertainty on i?^ 

The main sources of systematic uncertainties in the deter- 
mination of R^^-^ are related to the alignment accuracy of 
the PT system and to the determination of the f] value. 

The systematic uncertainty due to misalignment ef- 
fects was evaluated in different ways. A given offset Acj) — 
Aip+6(j) can be directly propagated in the algorithm which 
computes the charge ratio to evaluate ^ + SRfj,. 
The S(j) = 0.2 mrad uncertainty on the alignment ac- 
curacy obtained with magnets off (Sec. 13.51) corresponds 
to SRfi — 0.03. However a more powerful procedure was 
used to better estimate this systematics. We considered 
all muon tracks crossing both arms of each spectrome- 
ter, thus providing two independent deflection values Acjj 
per spectrometer for the same muon track. With perfect 
alignment and neglecting the energy loss the difference 
SAip = Acjjarmi — A(/)arm2 should be peaked at zero. The 
two distributions, one for each spectrometer, are shown in 
Fig. [5] together with a Gaussian fit to the central part of 
the distributions, where the effects of muon energy loss in 
the magnet iron are negligible. The two peaks are at 0.08 
mrad and -0.07 mrad respectively, ~2 standard deviations 
away from zero. A misalignment of 0.08 mrad produces 
an error on the charge ratio SRf^ ~ 0.015. We quote this 
number as the limiting alignment accuracy of each dou- 
blet with respect to the other. This number is conservative 
since it assumes that all four arms are affected indepen- 
dently from the same uncertainty. In reality only the the 
outer two doublets of each magnet contribute to this er- 
ror, since a given offset in the central doublet cancels the 
systematic uncertainty for Z\0armi and A(j)arm2- 

A further test which also incorporates local effects con- 
sists in comparing the values i?^ (j=l,...,4) in each magnet 
arm. The average difference from the mean value \R^^ — 
= 0.017 is within the statistical accuracy of each 
(5i?' = 0.03. 

Another consistency check exploited a small data sam- 
ple (~9 days of livetime) obtained after inverting the po- 
larity of the magnetic field. Running with inverted mag- 
netic polarity could in principle cancel the systematic error 
related to misalignment effects. The result is /j™''erted _ 
1.36 ± 0.04, corresponding to the unfolded value Jij^'"^^t<^d 
= 1.39 ± 0.04. Even if the statistical error is larger than 
the systematic error quoted above, the result is in good 
agreement with the value obtained with normal polarity. 

The charge- misidentification rj was previously estimated 
using Monte Carlo simulations. As already discussed the 
value is larger than what is expected from multiple scat- 
tering alone. The difference is ascribed to the inclusion of 
spurious effects, such as the production of secondary par- 
ticles near the muon trajectory, timing errors, and other 



second order effects not reproducible with the Monte Carlo 
program. Therefore we expect that the systematic uncer- 
tainty on rj is one-sided, being rjreai > Vmc- To estimate 
this difference rj was evaluated using experimental data 
for a subsample of events. We considered all muon tracks 
crossing both arms of each spectrometer, which provide 
two independent defiections Acj) of the same muon track. 
In this case, the probability that the two deflection an- 
gles have opposite sign is p — 277(1 — 77) and therefore 
7; = 1 — ^/l — 2p. This formula neglects the correlation 
between the two A(j) angles, since they are built using 
a common track (the one in between the two arms). The 
correct ri{p) relation was derived using a Monte Carlo sim- 
ulation applied to the experimental and simulated data. 
It was found for the case of doublets rjdata — 0.018 ±0.002 
and r]MC — 0.012±0.002. Considering doublets and mixed 
configuration together, we found rjdata = 0.026±0.002 and 
ryjvfc = 0.019 ± 0.002. The difference 6t] = 0.007 was used 
as an estimate of the systematic uncertainty on rj which 
corresponds to SR^ — 0.007. 

The final systematic error is taken as the quadratic 
sum of its contributions and it is assumed be the same for 
single and for multiple muon events: 

SR;-f i^yst.) =troll (13) 
5 i?^ as a function of and S^^ cos 9* 

The underground muon momentum was computed us- 
ing Eq. [51 The muon charge ratio as a function of p^ is 
shown in Fig. [BJ where the widths of the horizontal er- 
ror bars correspond approximately to the (average) muon 
momentum resolution. A linear fit 

RM^) - flo + ailogio[p^/(GeV/c)] (14) 

gives ao = 1.29 ±0.06 and ai = 0.05 ±0.03 with x^/dof = 
13.7/15. The data are also compatible with the hypothesis 
of a constant charge ratio, since the fit to a constant yields 
an = 1.379 ± 0.015 with x^/dof = 16.2/16 and therefore 
Ax^/dof — 2.47/1 (corresponding to ~1.6 sigma). 

The muon energy at the surface (f^) is directly related 
to the underground residual energy (E"^ — Pfi) and to the 
rock amount crossed by the muon to reach the detector 
level. In fact, the energy loss of high energy muons in the 
rock is usually expressed as 

HE 

- — ^ a{E) + I3{E)E (15) 

where h is the rock depth while the two energy-dependent 
parameters a and /3 are the contributions of the ionization 
energy loss and the radiative processes, respectively. Eq. 
[TSl can be integrated to obtain the approximate formula 

£^ = {E^ + a/l3)e'"' - a/l3 (16) 

which connects the surface and underground muon ener- 
gies. However, Eq. (THlis valid only on average. The "reso- 
lution" dSfj^ = EJ^'^'^ — is dominated by the statistical 
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V 


nun} 


Doublets 


13595 


9993 


1.360 ± 0.018 


0.0165 ± 0.0012 


1.375 ± 0.019 


Mixed 


8951 


6603 


1.355 ± 0.022 


0.0403 ± 0.0022 


1.393 ± 0.025 


Singlets 


2181 


1704 


1.28 ± 0.064 


0.064 ± 0.005 


1.33 ± 0.05 



Table 2: Final statistics for the underground muon charge ratio. Results are given separately for the three classes of 
events defined in the text. Errors are statistical only. 



(A) (gM)„„,„„.. H fraction jV,/jV„ J?-^ 

=1 3.35±0.09 (19.4±0.1) TeV 0.667±0.007 4.99±0.05 1.377±0.014 
>1 8.5±0.3 (77±1) TeV 0.352±0.012 2.09±0.07 1.23±0.06 

Table 3: Primary cosmic ray information for single and multiple muon events (see text). Reported numbers were 
obtained with MC2 and with the composition model fitted in ^D]. Only statistical errors are quoted. Systematic 
uncertainties related to the composition model dominate and can be inferred from the cited reference {S{A) ~ 1). In 
the last column the measured (and unfolded) charge ratios are given. 



Bin 


f ,1 cos 6* range 










SR";,"^ (stat) 


5i?r-'(sys) 




(GeV) 


(GeV) 










(%) 


1 


891 - 1259 


1197 


899 


678 


1.353 


0.074 


0.4 


2 


1259 - 1778 


1527 


14125 


10571 


1.373 


0.019 


0.5 


3 


1778 - 2512 


2071 


10345 


7613 


1.420 


0.025 


1.3 


4 


2512 - 3548 


2852 


3232 


2444 


1.409 


0.047 


4.3 


5 


3548 - 7079 


4329 


643 


548 


1.192 


0.079 


3.1 



Table 4: Main information for the five bins in £pCos0*. From left to right: the energy range and average value, the 
number of muons reconstructed with positive and negative charges, the unfolded charge ratio, the statistical and 
systematic errors. 



fluctuations due to the discrete processes described by the 
term /3 [21] ■ We evaluated f ^ with a full Monte Carlo sim- 
ulation to build the table £^ = f{h,p^). For this purpose 
the code MC2 was used since it contains a detailed de- 
scription of the muon flux at the surface and the muon 
transport in the Gran Sasso rock. The (/i,p^) plane was 
divided into 10x10 equally-spaced bins and in each bin 
the average {E^) value was computed. The binning was 
chosen coarse enough to have a large statistical sample 
in each bin without affecting the resolution d£^^ which 
is of the order of 0.15 in the logarithmic scale. The sur- 
face muon charge ratio was computed as a function of the 
variable (£^i) cos 9* binned according to the resolution. Fi- 
nally, the experimental values were corrected in each bin 
for the corresponding charge mis-identification and shown 
in Fig. [7] with black points (only single muon events were 
considered) . The present statistics does not allow to draw 
conclusions about the highest energy data point shown in 
the figure. 

Tab. Ogives some information for each of the five bins 
considered: the energy range and average value, the statis- 
tical sample, the unfolded charge ratio, the statistical and 
systematic errors. The latter were evaluated computing in 
each bin the two contributions discussed in Sec. 14.21 

In Fig. [7] are shown for comparison the data from other 
experiments for which we could recover information on the 
£fi cos 9* variable. For the low energy region we took data 
from Ref. [22| and Ref. [23] (we choose data points with 
uncertainties (5ii!^<0.05) while in the high energy region 
the data are from Ref. [S] and Ref. [TIT. For the latter, since 



the angular information were not provided in the paper, 
we plotted the i?^ integrated value in correspondence of 
the £fj, cos 9* value given in Ref. [^ . We also report a recent 
result from Ref. [24] where the vertical muon charge ratio 
is given in the range 1-3 TeV (average value 1.3 TeV). 

Finally, we fit our data to Eq. [S] using a procedure 
similar to what is described in Ref. [10] . We rewrite Eq. [5] 
in the form: 



1 + bTrC COS 9* /e^ I + bKC^ cos 9* /ck 

where Rktt = Znk/Zn-k and = 1-/^- = Zn^+/Zn-k 
(and similarly for kaons). f,^+ and fx+ were left free to 
vary while we fixed the kinematical parameters = 0.674, 
UK — 0.246, — 1.061, Bk = 1.126 and the fraction of 
kaons over pions in the atmosphere Rktv = 0.149 [5]. The 
fit of i?^ = /(j)^- takes into account data from [35] 
and [23] for the low energy region and data from this 
work at higher energies. The fit yields the values /7r+ 
= 0.5514±0.0014 and /a'+ 0.680±0.015 which corre- 
spond to a ratio R^r = Zj^j^+ /Z^^^- — 1.229±0.001 and 
Rk = Zfqj^+jZfqj^- = 2.12±0.03 for pions and kaons re- 
spectively. The result of the fit is shown in Fig. [7] as a 
continuous line. 

The contribution of the prompt muon component to 
i?^ was evaluated for three different charm production 
models: the phenomenological non-perturbative models 
RQPM and QGSM [2] and the semi-empirical model from 
Volkova et al. |25] . In [2] the prompt muon fiux and charge 
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Fig. 5: Two-arm test. Distributions of the difference of the deflection angles for tracks crossing both arms of one 
spectrometer: SMI (left) and SM2 (right). In each plot we show the fit of the central part of the distributions to a 
Gaussian function. 



ratios are parametrized as a function of the muon energy. 
The results of the fit extended to include the prompt con- 
tribution as predicted by these models, are shown in Fig. 
[71 The pion and kaon charge ratios obtained from the fit 
are unchanged within the statistical errors. 



6 Conclusions 

The atmospheric muon charge ratio i?^ — l^^i- 
measured using the spectrometers of the OPERA under- 
ground detector. We analyzed four months of data taken 
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Fig. 6: Measured charge ratio of underground muons as 
a function of the reconstructed muon momentum. Data 
points below ^5 GeV/c and above ^1000 GeV/c are sup- 
pressed by the cut on Z\0. A fit of the form i?p(p^) — 
oo -I- ailogio[p;j/(GeV/c)] is superimposed to the data. 



during the 2008 CNGS neutrino run. For single muons the 
Rfi value integrated over the underground muon spectrum 
is 

RT' ("m = 1) = 1-377 ± 0.014 (stat.)t°:°J^ (syst.) 

to be compared to R'^'^^inf, > 1) = 1.23 ± 0.06 for muon 
bundles. This difference of about ^2Ag supports the hy- 
pothesis of the decrease of the muon charge ratio with in- 
creasing primary mass. This is the first indication of such 
an effect which provides a further handle for the correct 
understanding and modelling of the secondary production 
in the atmosphere. 

The underground muon charge ratio is consistent with 
past measurements in a similar energy region. Data sug- 
gest a slight increase of i?^ with the underground muon 
momentum, although a fit to a constant charge ratio can- 
not be excluded. 

The dependence of i?^ on the vertical surface energy 
cos shows an increase in the region 1-3 TeV and it is 
compatible with a model which considers only the tt and 
K contributions to the muon charge ratio. A fit of the low 
energy data and our data with a simplified description of 
the atmospheric muon flux provides a value of the pion 
and kaon charge ratios Rt^ = Zj^^+ /Z^^^- ~ 1.229±0.001 
and Rk = Zf^K+ /Znk- — 2.12±0.03, respectively. The 
inclusion of the prompt muon component does not modify 
the fit results. It is however intriguing to observe that our 
measurement lies in the region where the charmed particle 
production may start to give an observable contribution to 
the muon charge ratio. A larger statistical sample or an 
experimental measurement with a new detector at very 
large depths could shed light on the region £^ cos 9* > 
10 TeV. The data collected by OPERA at the end of its 
scientific program will allow to improve the measurement 
in this energy region. 
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Fig. 7: Rf^ values measured by OPERA in bins of £f^ cos 9* (black points). Also plotted are the data in the low energy 
region from MINOS-ND ^ and L3+C t23j and in the high energy region from Utah ^ , MINOS [lU] and LVD [H] 
experiments. The result of the fit of OPERA and L3+C data to Eq. [17] is shown by the continuous line. The dashed, 
dotted and dash-dot lines are, respectively, the fit results with the inclusion of the RQPM, QGSM [2] and VFGS [IS] 
models for prompt muon production in the atmosphere. 
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A Unfolded charge ratio 

Let us call m*^ the number of muons with charge i recon- 
structed with charge j. The total number of true positive 
and negative muons is therefore: 

= m + m ^ 

On the other hand, the total number of reconstructed pos- 
itive and negative muons is: 

M+ = m++ + m"+ 
M" — m + 
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Let us define the charge-misidentification rj as: 



References 



+- _ 



m 



M- 



(18) 



Using a matrix notation, we can express the relationship 
between M and M as: 



M = HM 



where 



1 



■q- 



(19) 



(20) 



Inverting this relation, one has the number of "true" pos- 
itive and negative muons: 



M = H^M 



(21) 



where 



Hi 1 fl 

1 - ri+- - 7]-+ y - 



The two T] values r]^ and r/ are obtained from a Monte 
Carlo simulation. We found that, within the statistical 
accuracy of the simulation, 77 ^ = rj ^ = 77 as one would 
expect from a charge-symmetric detector. This simplifies 
the expressions which, in terms of the ratio R, becomes 



R 



M+ _ (1 - rj)M+ - r]M- 
M- ~ -TjA!+ + (1 - T])M- 



^■qR + (1 - 77) 



(23) 



where R = M+ / M' . 

If R is computed with the same Monte Carlo events 
used to evaluate rj, one would obtain the same "true" R 
value of the starting data sample. If R is computed with 
the experimental reconstructed data, then R is the un- 
folded experimental value in Eq. 1101 

The error 5R is obtained propagating the errors on R 
and 77 over Eq. 1231 
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a/(1 - 277)2(5i?)2 + (R2 _ 1)2(^^)2 

5R = ^ ^ -— — (24) 

[vR-{i-vW 



It may be pointed out that we did not use any regular- 
ization scheme in the unfolding, i.e. statistical fluctuations 
on R are not damped in Eq. [221 in order to prevent un- 
physical spikes in the unfolded R value. This is acceptable 
in our case since the collected statistics on M+ and 
is large enough. 



